"""
Comprehensive test suite for differential operators.

This module provides validation tests for all differential operators to ensure:
- Mathematical accuracy against analytical solutions
- Proper handling of boundary conditions
- Convergence rates for finite difference schemes
- Performance benchmarks
- Singularity handling in curved coordinates
- Conservation properties

Test categories:
1. Unit tests for individual operators
2. Accuracy tests against known analytical solutions  
3. Convergence rate validation
4. Boundary condition tests
5. Performance benchmarks
6. Integration tests with PDE solvers
"""

using Test
using LinearAlgebra
using Random

# Include all operator modules
include("common.jl")
include("cartesian.jl")
include("spherical_sym.jl") 
include("cylindrical_sym.jl")
include("polar_sym.jl")
include("boundary_conditions.jl")
include("optimization.jl")

# Export test functions
export test_all_operators, test_operator_accuracy, test_boundary_conditions
export test_convergence_rates, benchmark_operators, test_conservation
export test_singularity_handling

"""
    test_all_operators()

Run complete test suite for all operators.
"""
function test_all_operators()
    println("Running comprehensive operator tests...")
    
    @testset "Differential Operators Test Suite" begin
        test_cartesian_operators()
        test_spherical_operators()
        test_cylindrical_operators()
        test_polar_operators()
        test_boundary_conditions()
        test_optimization_features()
        test_conservation_properties()
    end
    
    println("All operator tests completed!")
end

"""
    test_cartesian_operators()

Test Cartesian coordinate operators.
"""
function test_cartesian_operators()
    @testset "Cartesian Operators" begin
        # Create test grids
        grid_1d = CartesianGrid((-1.0, 1.0), 50)
        grid_2d = CartesianGrid((-1.0, 1.0), (-1.0, 1.0), (50, 50))
        grid_3d = CartesianGrid((-1.0, 1.0), (-1.0, 1.0), (-1.0, 1.0), (20, 20, 20))
        
        test_laplacian_accuracy_1d(grid_1d)
        test_laplacian_accuracy_2d(grid_2d)
        test_gradient_accuracy(grid_2d)
        test_divergence_accuracy(grid_2d)
        test_curl_accuracy(grid_2d)
        test_biharmonic_operator(grid_2d)
        test_fft_operators(grid_2d)
    end
end

"""
    test_laplacian_accuracy_1d(grid)

Test 1D Laplacian accuracy against analytical solution.
"""
function test_laplacian_accuracy_1d(grid)
    @testset "1D Laplacian Accuracy" begin
        # Test function: u(x) = sin(πx)
        # Analytical Laplacian: ∇²u = -π²sin(πx)
        
        x_coords = coordinate_arrays(grid)[1]
        u_exact = sin.(π * x_coords)
        laplacian_exact = -π^2 * sin.(π * x_coords)
        
        # Create Laplacian operator
        laplace_op = make_laplace(grid)
        u_numerical = similar(u_exact)
        laplace_op(u_exact, u_numerical)
        
        # Compare interior points (exclude boundaries)
        interior = 2:(length(u_exact)-1)
        error = norm(u_numerical[interior] - laplacian_exact[interior], Inf)
        
        @test error < 0.1  # Should be much smaller with proper BC handling
        
        println("1D Laplacian max error: $error")
    end
end

"""
    test_laplacian_accuracy_2d(grid)

Test 2D Laplacian accuracy.
"""
function test_laplacian_accuracy_2d(grid)
    @testset "2D Laplacian Accuracy" begin
        # Test function: u(x,y) = sin(πx)sin(πy)
        # Analytical Laplacian: ∇²u = -2π²sin(πx)sin(πy)
        
        x_coords, y_coords = coordinate_arrays(grid)
        u_exact = [sin(π*x)*sin(π*y) for x in x_coords, y in y_coords]
        laplacian_exact = [-2π^2*sin(π*x)*sin(π*y) for x in x_coords, y in y_coords]
        
        laplace_op = make_laplace(grid)
        u_numerical = similar(u_exact)
        laplace_op(u_exact, u_numerical)
        
        # Compare interior points
        interior_x = 2:(size(u_exact,1)-1)
        interior_y = 2:(size(u_exact,2)-1)
        error = norm(u_numerical[interior_x, interior_y] - 
                    laplacian_exact[interior_x, interior_y], Inf)
        
        @test error < 0.5  # Adjust tolerance based on grid resolution
        
        println("2D Laplacian max error: $error")
    end
end

"""
    test_gradient_accuracy(grid)

Test gradient operator accuracy.
"""
function test_gradient_accuracy(grid)
    @testset "Gradient Accuracy" begin
        # Test function: u(x,y) = x²y + y²x
        # Analytical gradient: ∇u = (2xy + y², x² + 2yx)
        
        x_coords, y_coords = coordinate_arrays(grid)
        u_exact = [x^2*y + y^2*x for x in x_coords, y in y_coords]
        
        grad_x_exact = [2*x*y + y^2 for x in x_coords, y in y_coords]
        grad_y_exact = [x^2 + 2*y*x for x in x_coords, y in y_coords]
        
        # Test x-gradient
        grad_x_op = make_gradient(grid, 1)
        grad_x_numerical = similar(u_exact)
        grad_x_op(u_exact, grad_x_numerical)
        
        # Test y-gradient  
        grad_y_op = make_gradient(grid, 2)
        grad_y_numerical = similar(u_exact)
        grad_y_op(u_exact, grad_y_numerical)
        
        # Compare interior points
        interior_x = 2:(size(u_exact,1)-1)
        interior_y = 2:(size(u_exact,2)-1)
        
        error_x = norm(grad_x_numerical[interior_x, interior_y] - 
                      grad_x_exact[interior_x, interior_y], Inf)
        error_y = norm(grad_y_numerical[interior_x, interior_y] - 
                      grad_y_exact[interior_x, interior_y], Inf)
        
        @test error_x < 0.1
        @test error_y < 0.1
        
        println("Gradient errors - x: $error_x, y: $error_y")
    end
end

"""
    test_divergence_accuracy(grid)

Test divergence operator accuracy.
"""
function test_divergence_accuracy(grid)
    @testset "Divergence Accuracy" begin
        # Test vector field: v = (x²y, xy²)
        # Analytical divergence: ∇·v = 2xy + 2xy = 4xy
        
        x_coords, y_coords = coordinate_arrays(grid)
        nx, ny = size(grid)
        
        vector_field = zeros(nx, ny, 2)
        for (i, x) in enumerate(x_coords), (j, y) in enumerate(y_coords)
            vector_field[i, j, 1] = x^2 * y  # vx component
            vector_field[i, j, 2] = x * y^2  # vy component
        end
        
        divergence_exact = [4*x*y for x in x_coords, y in y_coords]
        
        div_op = make_divergence(grid)
        div_numerical = zeros(nx, ny)
        div_op(vector_field, div_numerical)
        
        # Compare interior points
        interior_x = 2:(nx-1)
        interior_y = 2:(ny-1)
        error = norm(div_numerical[interior_x, interior_y] - 
                    divergence_exact[interior_x, interior_y], Inf)
        
        @test error < 0.5
        
        println("Divergence max error: $error")
    end
end

"""
    test_curl_accuracy(grid)

Test curl operator accuracy.
"""
function test_curl_accuracy(grid)
    @testset "Curl Accuracy" begin
        # Test vector field: v = (-y, x) (rotation field)
        # Analytical curl: ∇×v = ∂vx/∂y - ∂vy/∂x = -1 - 1 = -2
        
        x_coords, y_coords = coordinate_arrays(grid)
        nx, ny = size(grid)
        
        vector_field = zeros(nx, ny, 2)
        for (i, x) in enumerate(x_coords), (j, y) in enumerate(y_coords)
            vector_field[i, j, 1] = -y  # vx component
            vector_field[i, j, 2] = x   # vy component
        end
        
        curl_exact = fill(-2.0, nx, ny)
        
        curl_op = make_curl(grid)
        curl_numerical = zeros(nx, ny)
        curl_op(vector_field, curl_numerical)
        
        # Compare interior points
        interior_x = 2:(nx-1)
        interior_y = 2:(ny-1)
        error = norm(curl_numerical[interior_x, interior_y] - 
                    curl_exact[interior_x, interior_y], Inf)
        
        @test error < 0.1
        
        println("Curl max error: $error")
    end
end

"""
    test_spherical_operators()

Test spherical coordinate operators.
"""
function test_spherical_operators()
    @testset "Spherical Operators" begin
        grid = SphericalSymGrid(50, (0.1, 2.0))
        
        test_spherical_laplacian(grid)
        test_spherical_singularity_handling(grid)
        test_spherical_conservation(grid)
    end
end

"""
    test_spherical_laplacian(grid)

Test spherical Laplacian accuracy.
"""
function test_spherical_laplacian(grid)
    @testset "Spherical Laplacian" begin
        # Test function: u(r) = r²
        # Analytical spherical Laplacian: ∇²u = d²u/dr² + (2/r)du/dr = 2 + 2*2 = 6
        
        r_coords = coordinate_arrays(grid)[1]
        u_exact = r_coords.^2
        laplacian_exact = fill(6.0, length(r_coords))
        
        laplace_op = make_laplace_spherical(grid)
        u_numerical = similar(u_exact)
        laplace_op(u_exact, u_numerical)
        
        # Compare interior points (avoid r=0 and boundary)
        interior = 3:(length(u_exact)-2)
        error = norm(u_numerical[interior] - laplacian_exact[interior], Inf)
        
        @test error < 0.1
        
        println("Spherical Laplacian max error: $error")
    end
end

"""
    test_spherical_singularity_handling(grid)

Test proper handling of r=0 singularity in spherical coordinates.
"""
function test_spherical_singularity_handling(grid)
    @testset "Spherical Singularity Handling" begin
        # Test function that's regular at r=0: u(r) = r²
        r_coords = coordinate_arrays(grid)[1]
        u_test = r_coords.^2
        
        # The Laplacian should give finite result at r=0
        laplace_op = make_laplace_spherical(grid)
        result = similar(u_test)
        laplace_op(u_test, result)
        
        @test isfinite(result[1])  # Should not be NaN or Inf at r=0
        @test abs(result[1] - 6.0) < 1.0  # Should be close to analytical value
        
        println("Value at r=0: $(result[1]) (expected: 6.0)")
    end
end

"""
    test_convergence_rates()

Test convergence rates of finite difference schemes.
"""
function test_convergence_rates()
    @testset "Convergence Rates" begin
        # Test second-order schemes should have O(h²) convergence
        resolutions = [25, 50, 100, 200]
        errors = Float64[]
        
        for n in resolutions
            grid = CartesianGrid((-1.0, 1.0), n)
            x_coords = coordinate_arrays(grid)[1]
            
            # Test function: u(x) = x⁴  
            # Exact Laplacian: u''(x) = 12x²
            u_test = x_coords.^4
            laplacian_exact = 12 * x_coords.^2
            
            laplace_op = make_laplace(grid)
            result = similar(u_test)
            laplace_op(u_test, result)
            
            # Compute error on interior points
            interior = 3:(length(result)-2)
            error = norm(result[interior] - laplacian_exact[interior], Inf)
            push!(errors, error)
        end
        
        # Check convergence rate
        for i in 2:length(errors)
            ratio = errors[i-1] / errors[i]
            expected_ratio = (resolutions[i] / resolutions[i-1])^2  # O(h²)
            
            # Allow some tolerance in convergence rate
            @test ratio > expected_ratio * 0.5
            @test ratio < expected_ratio * 2.0
            
            println("Resolution $((resolutions[i-1], resolutions[i])): ratio = $ratio, expected ≈ $expected_ratio")
        end
    end
end

"""
    test_boundary_conditions()

Test boundary condition implementations.
"""
function test_boundary_conditions()
    @testset "Boundary Conditions" begin
        grid = CartesianGrid((0.0, 1.0), 50)
        
        test_dirichlet_bc(grid)
        test_neumann_bc(grid)
        test_robin_bc(grid)
        test_periodic_bc()
    end
end

"""
    test_dirichlet_bc(grid)

Test Dirichlet boundary condition implementation.
"""
function test_dirichlet_bc(grid)
    @testset "Dirichlet BC" begin
        # Create boundary condition set
        bc_left = DirichletBC(1.0)
        bc_right = DirichletBC(0.0)
        bcs = [bc_left, bc_right]
        
        # Create Laplacian matrix
        matrix, rhs = make_laplace_matrix_cartesian(grid, bcs)
        
        # Check that boundary rows are set correctly
        n = size(matrix, 1)
        
        # First row should be [1, 0, 0, ...]
        @test matrix[1, 1] ≈ 1.0
        @test sum(abs.(matrix[1, 2:end])) < 1e-12
        @test rhs[1] ≈ 1.0
        
        # Last row should be [..., 0, 0, 1]  
        @test matrix[n, n] ≈ 1.0
        @test sum(abs.(matrix[n, 1:end-1])) < 1e-12
        @test rhs[n] ≈ 0.0
        
        println("Dirichlet BC implementation verified")
    end
end

"""
    test_conservation_properties()

Test conservation properties of operators.
"""
function test_conservation_properties()
    @testset "Conservation Properties" begin
        # Test that conservative operators preserve mass
        grid = SphericalSymGrid(50, (0.1, 2.0))
        
        # Create conservative Laplacian
        laplace_conservative = make_conservative_laplace_spherical(grid)
        
        # Test with constant function (should give zero)
        u_constant = ones(grid.num_r)
        result = similar(u_constant)
        laplace_conservative(u_constant, result)
        
        # Sum should be close to zero (conservation)
        total_result = sum(result)
        @test abs(total_result) < 1e-10
        
        println("Conservation test: sum of Laplacian of constant = $total_result")
    end
end

"""
    benchmark_operators()

Benchmark operator performance.
"""
function benchmark_operators()
    @testset "Performance Benchmarks" begin
        println("\nOperator Performance Benchmarks:")
        println("=" * 50)
        
        # Test different grid sizes
        sizes_1d = [100, 500, 1000, 5000]
        sizes_2d = [50, 100, 200]
        
        benchmark_1d_operators(sizes_1d)
        benchmark_2d_operators(sizes_2d)
    end
end

"""
    benchmark_1d_operators(sizes)

Benchmark 1D operators.
"""
function benchmark_1d_operators(sizes)
    println("\n1D Operator Benchmarks:")
    
    for n in sizes
        grid = CartesianGrid((-1.0, 1.0), n)
        test_data = randn(n)
        result = similar(test_data)
        
        # Benchmark Laplacian
        laplace_op = make_laplace(grid)
        time = @elapsed for _ in 1:100
            laplace_op(test_data, result)
        end
        
        throughput = 100 * n / time  # operations per second
        println("  n=$n: $(round(time*10, digits=2))ms ($(round(throughput/1e6, digits=2))M pts/sec)")
    end
end

"""
    benchmark_2d_operators(sizes)

Benchmark 2D operators.
"""
function benchmark_2d_operators(sizes)
    println("\n2D Operator Benchmarks:")
    
    for n in sizes
        grid = CartesianGrid((-1.0, 1.0), (-1.0, 1.0), (n, n))
        test_data = randn(n, n)
        result = similar(test_data)
        
        # Benchmark Laplacian
        laplace_op = make_laplace(grid)
        time = @elapsed for _ in 1:10
            laplace_op(test_data, result)
        end
        
        total_points = n * n
        throughput = 10 * total_points / time
        println("  $(n)×$(n): $(round(time*100, digits=2))ms ($(round(throughput/1e6, digits=2))M pts/sec)")
    end
end

"""
    test_optimization_features()

Test optimization and caching features.
"""
function test_optimization_features()
    @testset "Optimization Features" begin
        grid = CartesianGrid((-1.0, 1.0), (-1.0, 1.0), (50, 50))
        bcs = [DirichletBC(0.0), DirichletBC(0.0), DirichletBC(0.0), DirichletBC(0.0)]
        
        # Test operator caching
        cached_op1 = cache_operator(() -> make_laplace(grid), grid, bcs, :test_laplacian)
        cached_op2 = get_cached_operator(grid, bcs, :test_laplacian)
        
        @test cached_op2 !== nothing
        @test cached_op1.cache_key == cached_op2.cache_key
        
        # Test memory pool
        temp_array = get_temp_array(Float64, (50, 50))
        @test size(temp_array) == (50, 50)
        
        return_temp_array!(temp_array)
        
        # Get cache statistics
        stats = get_cache_statistics()
        @test stats.cached_operators >= 1
        
        println("Optimization features verified")
    end
end

"""
    run_integration_tests()

Run integration tests with actual PDE problems.
"""
function run_integration_tests()
    @testset "Integration Tests" begin
        # Test with diffusion equation: ∂u/∂t = D∇²u
        test_diffusion_integration()
        
        # Test with Poisson equation: ∇²u = f
        test_poisson_integration()
    end
end

"""
    test_diffusion_integration()

Test operators in context of diffusion equation.
"""
function test_diffusion_integration()
    @testset "Diffusion Integration" begin
        grid = CartesianGrid((0.0, 1.0), (0.0, 1.0), (50, 50))
        
        # Initial condition: Gaussian
        x_coords, y_coords = coordinate_arrays(grid)
        initial_condition = [exp(-10*((x-0.5)^2 + (y-0.5)^2)) 
                           for x in x_coords, y in y_coords]
        
        # Diffusion coefficient
        D = 0.1
        dt = 0.001
        
        # Create Laplacian operator
        laplace_op = make_laplace(grid)
        
        # Simple forward Euler step
        current_state = copy(initial_condition)
        laplacian_result = similar(current_state)
        
        for _ in 1:10  # Few time steps
            laplace_op(current_state, laplacian_result)
            current_state .+= dt * D * laplacian_result
        end
        
        # Check that solution is still positive and bounded
        @test all(current_state .>= 0)
        @test maximum(current_state) <= maximum(initial_condition)
        
        println("Diffusion integration test passed")
    end
end